D-A038 775 NAVAL POSTGRADUATE SCHOOL MONTEREY CALIF F/6 18/11 
FINITE ELEMENT SOLUTION OF A THREE=DIMENSIONAL NONLINEAR REACTO=--ETC(U) 
DEC 76 E C BERMUDES 


UNCLASSIFIED 


i 1O «ms We 
5_ 


i a2 
gy 6.” 


Wee 


22 Wee Mes 


NAVAL POSTGRADUATE SCHOOL 


Monterey, California 


To) 
i~ 
~- 
oe) 
me) 
© 
=< 
Q 
<x 


THESIS 


FINITE ELEMENT SOLUTION OF A THREE- 
DIMENSIONAL NONLINEAR REACTOR 
DYNAMICS PROBLEM WITH FEEDBACK 
by 
Eulogio Conception Bermudes 


December 1976 


Thesis Advisor: D. Salinas 
Thesis Advisor: D. H. Nguyen 


Approved for public release; distribution unlimited. D D C 


APR 29 197 


2 | ee 
DDC FILE COPY 


f 
, 
, 


oo 
SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered) 


REPORT DOCUMENTATION PAGE 


ze; Finite Element Solution of a [hree- 
.4 = Dimensional Nonlinear Reactor 
Dynamics Problem with "Feedback. 


I ri 7. AUTHOR(e) 


Eulogio Conception /Bermudes 
o 9. PERFORMING ORGANIZATION NAME AND TGnREaE 
x Naval Postgraduate School 
Monterey, California 93940 


¥ . CONTROLLING OFFICE NAME ANO ACORESS 
Naval Postgraduate School 
Monterey, California 93940 
. MONITORING AGENCY NAmE &@ ACORESS(/{ different trom Controlling Office) 1S. SECURITY CLASS. (of thie report) 


Naval Postgraduate School Sa Unclassified 
Monterey, California 93940” Gf 


YY 


ny one 


DISTRIBUTION STATEMENT (of thie Report) 


Approved for public release; distribution unlimited. 


. OISTRIBUTION STATEMENT (ef the abstract entered in Block 20, |{ different frem Report) | SESSION tor 
Aveke 


i 


. KEY WORDS (Continue on reveree side if neceeecary and identity by black number) 


. ABSTRACT (Continue on reveree side if neceecary and identity by bleck mumber) 

‘This work examines the three-dimensional dynamic response 
| of a nonlinear fast reactor with temperature-dependent feed- 
back and delayed neutrons when subjected to uniform and local 
| disturbances. The finite element method was employed to re- 
duce the partial differential reactor equation to a system of 
ordinary differential equations which can be numerically in- 
tegrated. A program for the numerical solution of large _ + 


| 
| 


a | ion 7s 1473 corTion oF | nov 6818 OMsoLere 


S/W 102-014-6601 | 
(Page 1) A tight ute [SECURITY CLASSIFICATION OF THis PAGE (When Dete Entered) 2, 


SNE 
Se CuMmry CLASSIFICATION OF THIS PAGE/(When Dota Entered. 


~X sparse systems of stiff differential equations developed by 
Franke and based on Gear's method solved the reduced neutron 
dynamics equation. Although a study of convergence by refin- 
ing element mesh.sizes was not carried out, the crude fi- 
nite element mesh.utilized yielded the correct trend of 
neutron dynamic behavior. 


DD Form, 1473 


1 Jan 72 
S/N 102-014-6601 2 SECURITY CLASSIFICATION OF THIS PAGECMREn Date Entered) 


FINITE ELEMENT SOLUTION OF A THREE- 
DIMENSIONAL NONLINEAR REACTOR 
DYNAMICS PROBLEM WITH FEEDBACK 


by 


Eulogio Conception Bermudes 
Lieutenant, United States Navy 
B.S., United States Naval Academy, 1970 


Submitted in partial fulfillment of the 
requirements for the degrees of 
MASTER OF SCIENCE IN MECHANICAL ENGINEERING 
and 


MECHANICAL ENGINEER 


from the 


NAVAL POSTGRADUATE SCHOOL 
December 1976 


Author Elegie Cmaxfliir, Bormualia— 


Approved by: ‘2 a 
1S visor 


Thesis Advisor 


had SAanhe, 


econd Reader 


REE OS. 


Chairman, Department o echanical Engineering 


U Jpop — 
Dean of Science and Engineering 
3 


TRA TEN EN REIN IIR 


ABSTRACT 


This work examines the three-dimensional dynamic response 
of a nonlinear fast reactor with temperature-dependent feed- 
back and delayed neutrons when subjected to uniform and local 
disturbances. The finite element method was employed to re- 


duce the partial differential reactor equation to a system of 


ordinary differential equations which can be numerically in- 


tegrated. A program for the numerical solution of large 
sparse systems of stiff differential equations developed by 
Franke and based on Gear's method solved the reduced neutron 
dynamics equation. Although a study of convergence by refin- 
ing element mesh sizes was not carried out, the crude fi- 
nite element mesh utilized yielded the correct trend of neu- 


tron dynamic behavior. 
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I. INTRODUCTION 


The nuclear reactor with delayed neutrons and tempera- 
ture-dependent feedback is a nonlinear system, whose response 
to uniform and local disturbances differs greatly from that 
of a linear reactor. The prompt temperature feedback model 
and one average group of delayed neutrons were incorporated 
in this work, Practically all neutron dynamics analysis to- 
day deals with two-dimensional geometry, implying a symmetry 
in neutron dynamics behavior during transient [1,2]. This 
Symmetry assumption, however, could be unrealistic in the 
safety analysis of nuclear reactors. A unique feature of 
this work is the consideration of the three-dimensional depen- 
dence of the neutron flux. No symmetry assumptions were im- 
posed on the problem. This work investigated neutron dynamics 
under uniform and local disturbances in the core. A uniform 
initial flux throughout the interior of the reactor was im- 
posed. The finite element method (FEM) was employed to solve 
this nonlinear initial-boundary value problem. The FEM is 
quite effective in handling discontinuous forcing functions 


thereby making it particularly suited for examining the ef- 


fects of localized perturbations and space-dependent feedbacks. 


Vent f a 


II. THE NUCLEAR REACTOR WITH TEMPERATURE 
DEPENDENT FEEDBACK AND DELAYED NEUTRONS 


A. THE PHYSICAL SYSTEM 

The system under consideration is a fast reactor of cy- 
lindrical geometry that is composed of two different regions. 
The reactor core or region I is cylinderical in shape and is 
fueled by U-235. Region II or the reflector region complete- 
ly surrounds the core and is composed of U-238. Both regions 
were assumed to be homogeneous. Table I lists the physical 
Properties and geometry for each region. A schematic of the 
fast reactor geomenry is shown in figure 1. 

Temperature and delayed neutron effects were taken into 
account. Also, a one-velocity or one-group model was assumed, 
thereby making the velocity independent of spatial or tem- 
poral effects. The delayed neutrons were considered only in 
region I since region II was assumed to be a non-multiplying 
medium. The temperature effects were also assumed to be only 
in the core region. 


In general form, the one-velocity neutron diffusion equa- 


tion is [3] 

aN(r,t) = 2 at iy. r 

+ = vDV°N - E.vN + S(r,t) (1) 
where r= (x,y,z) 


D = neutron diffusion coefficient 


=, = macroscopic absorption cross-section 
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number of neutrons in a volume element dV at 
a point r at time t 


ll 


TABLE I 


Physical Constants 


| SYMBOL DEFINITION VALUE 
Si radius of Region I 60 om 
Riz total reactor radius 30 cm 
Hy height of Region I 160 cm 
Hrr total reactor height 220 cm 
V neutron speed 4. 8x10’ cm/sec 
D: core neutron diffusion coefficient 0.913 cm 
Dry reflector neutron diffusion 1.200 cm © 
coefficient 
z core neutron absorption cross 0.01401 om 
=e section 
Ly reflector neutron absorption 0.0040 om + 
ai cross section 
1 v number of neutrons per fission 2.54 
critical fission cross section 0.005736 on 
8 delayed neutron fraction; 0.00642 
dollar reactivity 
\ - 
€ fission energy 7.652x10 
cal/ fission 
is Rc) modified convection heat 0.0632 3 & 
4 transfer coefficient cal/(cm"sec ~C) 
F a temperature coefficient -0.004/°C 
x abundance-weighted mean decay 0.4350 sec™2 


constant 
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a. top view 
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I - core b. front view 
II - reflector 


_ Figure 1. Schematic of cylindrical reactor. 
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vOV*NdV = number of neutrons diffusing into dV per 
unit time at time t 


zavN dV = number of neutrons absorbed in dV per 
unit time at time t 


S(r,t)dV = number of neutrons produced in dV per 
unit time at time t 
The neutron number density is related to the neutron flux by 


the expression [4] 
o(r,t) = vN(r,t) (2) 


where (r,t) is the flux at time t . The neutron diffu- 


sion equation in terms of the flux is depicted by 


} ao(r,t) = 2 2 y 
— oo = DV7g - £.6 + S(¥,t) (3) 


B. PROMPT AND DELAYED NEUTRONS 

The source or production term in equation (3) is composed 
of the contributions of the prompt and delayed neutrons. The 
majority of the fission neutrons are prompt neutrons that 
appear almost instantaneously (within 107? second) on fission. 
Assuming a fast neutron non-leakage probability of unity, the 


prompt neutron source is described by 


Sy (r,t) = (1-8)K (r,t) o(r,t) (4) 
where 
kK (t,t) = infinite medium multiplication factor 
8 = total fraction of delayed neutrons 


The ipecbien of delayed neutrons is very small (note 
that 8 = 0.00642 from Table I). However, they have a very 
significant effect on the reactivity because their mean life- 
times are long. Without delayed neutrons, reactor control 
would not be possible. These delayed neutrons are born in 
the decay by neutron emission of nuclei produced following 
the g-decay of certain fission fragments. For example, the 


87 86 


B-decay of the fission fragment Br leads to Kr 


87 


plus a 
neutron. Nuclei such as Br whose production in fission 
eventually leads to the emission of a delayed neutron are 
called delayed neutron precursors [4]. 

There are six main groups of delayed neutrons. Each 


group is classified according to its decay constant. The de- 


layed neutron source term is portrayed by 


where 


decay constant of the a group 


density of the yeh precursor 
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-— 

“St 

w 

ct » 
— 
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Assuming that the fission fragments do not migrate appre- 
ciable distances and assuming a non-circulating fuel reactor 
[3], the precursor density is delineated by 


aC. (r,t) 


aa = BK t,o - A5C, (6) 


th 


where 8. is the fraction of delayed neutrons of the i 


group. The solution to the precursor equation is in terms 


of a time integral expressed by 
t 
eth ty > ota f ey ltt) (Fe ty e(F t)dt! (7) 


fe) 
Inserting equation (7) in equation (5) yields the delayed 


neutron production term as 


k er 
Byte ff ev Alt (CF tyo(F,t)dt! 


- 6 
Sy(r.t) "2a 
i 0 (8) 


] 
For convenience the six main groups of delayed neutrons 
were considered as one group. This was accomplished by us- 


ing the abundance-weighted mean decay constant defined by 


6 
i=] 


This is a reasonable approximation since many of the impor- 
tant phenomena of nuclear reactor dynamics can be character- 
ized satisfactorily by combining all the emitters or precur- 
sors into one, two, or, at most, three effective groups [3]. 
Replacing A; with the abundance-weighted mean decay con- 


stant showed the delayed neutron source term to be 


t 
Sy(Fst) = eiz, f eA(t-t kK CR tyo(ryt)dt! 
© 
DOPPLER EFFECT AND TEMPERATURE FEEDBACK MODEL 
The Doppler effect in fast reactors is due to the tempera- 
ture broadening of many closely spaced high-energy resonances 


in both the fission and parasitic-absorption cross-sections. 
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These resonances basically mean that as the temperature in- 
creases, the number of neutrons that are absorbed increases 
and, similarly, the number of fission events increases. 
These nonproductive and productive processes compete in a 
complicated manner, and the net effect may be either an in- 
crease or decrease in reactivity [3]. For most fast reac- 
tors, the effect is negative, and it will be so assumed in 
this work. 

The reactivity change with respect to the fuel tempera- 


ture change is modeled by [2] 
+ crime) (11) 


where a, b, and c are parameters determined from experimen- 
tal or neutronic calculations, m is an integer, and T is 
the current fuel temperature. For most fast reactor cores 
(ceramic-fueled cores), TdK,/dT is very close to being con- 
stant over a wide range of temperatures. Therefore, it is 


assumed that a and c are zero and the Doppler coefficient 


is expressed by 


The reactivity model is then expressed by 


K_ (Ft) = Ko + btn (7) | (13) 
where 
Ko = multiplication factor at steady-state 
vy = number of neutrons produced per fission 
uP = original fuel temperatures 
b = Doppler constant 


The definition of K® is 


Vie 
° - 
i. ee (14) 
a 
where 
* 
Le = critical fission cross-section 
. = macroscopic absorption cross-section of the 


core 
The rise in fuel temperature is depicted by 


Sirah) * TlRsth = 7,08) (15) 


This is also described by the integral [5] 


t 
e(r,t) = f(t-t') w(r,t)dt' (16) 
J 


where f(t-t') is the feedback kernel and is dependent upon 


the type of temperature feedback model used. w(r,t) is the 


rise in neutron flux above the steady state and is expressed 


by 
v(r,t) = o(r,t) - 5 (4) (17) 
where d5(r) is the neutron flux at steady state. 


There are three types of temperature feedback models that 
can be considered here. The first is known as "Newton's 
feedback" which determines the reactor temperature by Newton's 
law of cooling. The second temperature feedback model is 
called the adiabatic feedback model. This represents the 
temperature for the loss of coolant case. The third model is 
the prompt temperature feedback model in which the fuel tem- 
perature follows the behavior of the neutron flux without de- 
lay [5,6]. The prompt feedback model was employed in this 
analysis. The feedback kernel for this temperature feedback 


model is 
f(t-t') = & 5(t-t') (18) 


where K is an energy production operator with units of °C 
per unit flux and y , a dimensionless quantity, is related 
to the mean time for heat transfer to coolant. 

Inserting equation (18) in equation (16) and performing 


the integration produced a rise in temperature of 


a(Ft) = © y(F,t) (19) 


STEM 


gies 


From equation (19) the temperature of the fuel is 
T(#,t) = 2 v(Ft) + T)(F) (20) 


The ratio of the current temperature to the steady-state tem- 


perature is therefore 


a ae a (21) 
ee ae 


In this work un was considered to be constant throughout 
the core. Incorporating equation (21) into equation (13) 


gave the reactivity model of 
K.(F,t) = ko + b on [—— v(F,t) + 1] (22) 
oo 0 ae 


D. FIELD EQUATIONS : 

Before establishing the field equations it is desirable 
to express equation (3) in terms of the rise in flux. Using 
equation (17) in equation (3) and grouping terms yielded the 
following diffusion equation: 


t 
1 ay = 2. J Ys J -X(t-t') ' 
v oe 7 COV y-tiy +(1-8)K 2.0 + BAL, e K vdt'] 
0) 


t 
2 - To -A(t-t') 1 
+ [DV*o.-2.o5+(1-B)K2,,+ BAL, f e Kot! ] 
(23) 


The second bracketed term in equation (23) is identically 


equal to zero since it is the steady state portion of the 


— 


diffusion equation. The rise in neutron flux above its 


steady state value is therefore expressed, for the core, by 


t 2H - 2 3 
a DV*y - Toy + (1-8) K24¥ 


t - ' 
+ Bi, J en Alt") gy dt! (24) 


fe) 


and, for the reflector region, by 

1 ov 2 pv2y - 

ae DV7p z,Y (25) 
Inserting the reactivity model into equation (24) yielded 


1 3v . py2y - ig ° 
: Dv?y - Ey + (1-8) 2AK2y 


ot 
v ; X ( ) 

ie K Tyo -i t=t" fn 

+(1-6)bz, taal + 1]} p+ ems f e ydt 
t y 1 
+ gizb if en A(t-t Tan + 1)] wdt'} (26) 
0 
(0) 


For the reflector, the last four terms of equation (26) are 
zero. The effects of the temperature on the delayed neutrons 
were neglected in this work. The field equations can now be ; 


expressed, for the core, as 


ie 
oy _ 2 - "i r 
st vDV7p + [vz v(1-8) vi gly ie 
fre Ky be 
+ [-(1 B)vz.b] Can (TT + 1)]y 3 
. © oX(t-t"} i 

+ [-Brviev] Cf e vdt'] = 0 (27) 
0 ee, 
3 
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and, for the reflector, as 
ow _ 2 m 
Ft vDV“wv + vay 0 (28) 


The non-linear terms of the core field equation will be 
linearized accordingly. In more compact form, equations (27) 


and (28) became 


t y 1 
ae . ctvtp + ozo + chfen(he eryiyteste 78h te? pace ten 
st YT, : 
(29) 


and 


_ - clv2y + c2y = 0 (30) 


where the meanings of the coefficients cl, c2, etc., are 
obvious for the core and reflector. 

Equations (29) and (30) were subjected to the following 
conditions: 

boundary condition: 


Up(F,.t) = 0 (31) 


where re are coordinates of points on the outer surface of 
the reflector and the subscript R refers to the reflector 
continuity of flux: 


va(r,t) = ¥.(Fp>t) (32) 


where ry are coordinates of points on the core-reflector 


interface and the subscript c refers to the core. 


III. FINITE ELEMENT 


A. INTRODUCTION 

The application of the finite elements of nonlinear con- 
tinua has been mostly in the field of solid mechanics. 

Prior work [1] using the finite elements of nonlinear con- 
tinua on a nonlinear reactor dynamics problem has been suc- 
cessful. The FEM in this work was utilized to reduce a non- 
linear partial differential equation of the nuclear reactor 
to a system of nonlinear ordinary differential equations in 
time. The time integration was accomplished by using a com- 
puter program for the numerical solution of stiff differen- 
tial equations developed by Franke [7]. In order to minimize 
computer storage requirements, an optimum compacting scheme 
(OCS), described by Ref. 8, was adopted. 

The finite element models of operator equations are gen- 
erally classified into three categories: a) the variational 
finite element models such as the Ritz method, b) the weigh- 
ted residuals method such as the method of Galerkin, and c) 
the direct finite element models which are not based on 
functional minimization. From experience in structural 
mechanics, the most effective method for generating accept- 
able finite element models of nonlinear equations is the “ 
Galerkin method [1]. This work adopted the method of Galer- 


kin in a finite element approximation over the spatial 


domain of the field equations. 


B. THE METHOD OF GALERKIN 


The Galerkin method is a special case of the method of 
weighted residuals. It involves a rational choice of 
weighting function that is consistent with the type of fi- 
nite element approximation considered. Indeed, the weight- 
ing functions chosen are the basis or shape functions em- 
ployed in the finite element approximation. There are two 
favorable characteristics of the Galerkin method which 
makes it attractive. The first attribute is its amenability 
to integration by parts. This supplied the freedom of using 
a lower order finite element than might be otherwise possible. 
The second favorable characteristic of the Galerkin method 
is that the symmetric operators in the field equations trans- 
form into symmetric matrix operators. Both these attributes 
are attractive for computational purposes. 


Consider the initial-boundary-value problem 


seedy - f(*,t) (33) 


where P contains the nonlinear operators. According to 
the spatial finite element discretization, the solution of 
equation (33) is in the form of an union of an N-term approx- 


imation given by 


where N is the number of system degrees of freedom (i.e., 
number of coordinates), and G(r) are the system or global 


basis functions which span the space of the approximate 


24 


solution w(r,t) [1]. The Einstein summation is used. The 
global basic functions are “pyramid functions", each of 
which has a prescribed functional description over a sub- 
domain of the system and is zero elsewhere [9]. The unknown 
coordinate functions y(t) are the time-dependent magni- 
tudes of the approximated flux wp(r,t) and/or its deriva- 
tives at discrete nodal points [10]. 


The residual function, R(r,t) , is defined such that it 
is identical to zero when (r,t) is equal to the exact solu- 
tion. The residual function is expressed by 


R(r,t) ee a (35) 


The Galerkin orthogonality condition (using the basis func- 
tions as weight functions), when applied to the residual func- 


tion, requires that 


R(r.t)dvol aU 4 ROT ie cig gh (36) 


Vol 


From the field equations, the residual function for the core 


TS i 
= bye OY rs i ky 4, 
R(r,t) st clV*y + c2w + alias kt 1) Ju 

t - | 

+ c5[ if en A(t-t ) vdt'] (37) if 

0 Ee 

and for the reflector 3 

R(7,t) = 24 - clv2h + c2y (38) 4 


ro 


Using equation (34) and applying Galerkin's orthogonality 


condition produced a core equation of 


Joetoysy(e)-v26,Ce1-¥1, + Gy(c2-y) | 


Vol 
Lan (7 w¥e + VW)I Gy(c4ey), 
t - 
tf ea A(t-t )(c5+y), dt']G,} = 0 (39) 
ie) 
where ee ey eee 


The reflector has a similar equation without the nonlinearity 


and is expressed by 


i 2 e; ee; } = 
feste,gcts - Vv Gj(cl v) 4 + G y(c2 v) 3) dVol 0 (40) 
Vol 


C. THE ELEMENT 

A three-dimensional quadratic isoparametric element was 
employed in this work. The parent element is a triangular 
prism or solid wedge with straight sides. The element shape 


functions are expressed in terms of area coordinates in the 


plane of the triangle and by an isoparametric coordinate 


along the prism axis. Figure 2 shows the parent element. 
This element was chosen because of the ease with which it 
fits the cylindrical structure when it is transformed into 

a curved element. This type of element has been used before 
as filler elements [11]. 

The area coordinates are defined by area ratios. Con- 
sider the triangle shown in figure 3. An arbitrary point P 
within the triangle defines three subareas designated by Als 
Ags and Ax. The ratio of each of the subareas to the total 
area is known as an area coordinate. In equation form, the 


area coordinates are 


L, = A,/A (41a) 

Ly = Ag/A (41b) 

Lz = A3/A (41¢c) 
where A _ is total area of the triangle. Ly» Los and E 


are the natural coordinates for a triangle. The requirement 
that the sum of the subareas be equal to the total area is 


obviously satisfied by the identity 
i. @ Eee ee (42) 


In. the plane of the triangle, only two of the area coordi- 
nates are independent. 

The isoparametric coordinates are best visualized by con- 
sidering the rectangular prism shown in figure 4. Isopara- 
metric coordinates are normalized coordinates such that 


their values on the faces of the rectangle are + 1. The Enc 


axes are in general not orthogonal. They are orthogonal 
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Figure 2. Quadratic triangular prism parent element 
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Figure 3. Definition of area coordinates 


Figure 4. Isoparametric coordinates 


only in the special case of a rectangular prism element 
[12]. The element basis functions use the ¢t coordinate 
in the prism axis and the Ly> Los and L3 coordinates in 
the plane of the triangle. 

The parent element, as shown in figure 2, has 15 local 
nodal points around its periphery. As such, there are 15 


basis functions which are given below [11]: 


Corner nodes: (nodes 1, 3, 5, 10, 12, 14) 
Ly(2Ly-1) (142) - 
Lo(2L,-1) (1+) 
L4(2L,-1)(1+2) 
Ly (2Ly-1) (1-2) 
Lo (2b 9-1) (1-5) 
L3(2L,-1) (1-2) 


Ny = 


Na = 
Ne = 
10° 
a 


Ny4* 


N 


Ni Ml ies tel SN SN 


Midside nodes of rectangles: (nodes 7, 8, 9) 
L, (1-27) 
= Lo(1-5*) 
= L4(1-¢*) 
Midside nodes of triangles: (nodes 2, 4, 6, 


= 2LjL,(1+s) 


= 2L,l,(1+c) 
= 2L,L, (1+) 


= 2L,L,(1-s) 
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The coordinates of each local node, in terms of Li» Los L3> 
and t are listed in table II. These element basis functions 
< N > define the geometry of the element. Note that they 


satisfy the relationship 


1 , at node i 


O , at other nodes (43) 


On the element level, the variation of the unknown func- 


tion wp is approximated by 
pe =< Nn '> {yp }® (44) 


where < N'> = row vector of element shape functions 


{wp }© = column vector of time-dependent nodal 


values of ¥° 


To satisfy continuity requirements, the shape functions 
< N'> have to be such that the continuity of the unknown 
function wp is preserved in the parent coordinates [11]. 
The shape functions < N > which characterize the element 
geometry and the shape functions < N'> which describe the 
unknown function do not necessarily have to be the same. 
There is no requirement that the nodal values be associated 
with the same nodes which were used to define the element 
geometry, though in practice it is often the case. Consider 
for example, the illustrations in figure 5, If the nodes 


defining the element geometry and the nodes defining the 


TABLE II 


Coordinates of Local Nodal Points 


Local Node 
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(a) Isoparametric 


(b) Super~parametric 


(c) Sub-parametric 
@ - geometry nodes 


O - variable function nodes “2 


Figure 5. Element classification 


unknown function are identical, the element is known as an 
isoparametric element. This means that the shape functions 
describing the geometry and the shape functions describing 


the unknown function w are equal, or 
<t'> «=< > (45) 


If there are more nodes defining the geometry than nodes de- 
fining the variable function, the element is called a super- 
parametric element. Using more nodal points to define the 
unknown function than to describe the geometry of the element 
results in a subparametric element [11]. This work utilized 


the isoparametric element classification. 


D. DIVISION OF THE SYSTEM INTO ELEMENTS 

In three-dimensional space, the division of the sys- 
tem into discrete elements is difficult to visualize. It is 
virtually impossible to show every nodal point of the system 
in one schematic. To present a clear view of the discretized 
domain, a “layer” approach was adopted. The first finite 
element grid or mesh employed here consisted of 128 elements. 
Under this grid (mesh I), the reactor was divided into four 
layers as shown in figure 6. Each layer was composed of 32 
elements. The first and fourth layers were each 30 cm in 
height and each contained entirely reflector elements. The 
second and third layers were each 80 cm in height and together 
they encompassed the entire core plus the remaining refiector 


elements. Each layer, in turn, was partitioned into three 


horizontal (xy) planes. The top plane included all the global 


— Z=110 


— 2=80 


Layer 2 


— 2=-80 


Layer 4 


—Z=-110 
I - core II - reflector 


Figure 6. Layers of mesh I 


or system nodes corresponding to local or element nodes 1 


through 6. The middle plane contained all the system nodes 
corresponding to the local nodes 7, 8, and 9. System nodes 
corresponding to element nodal points 10 through 15 comprised 
the bottom plane. To fix ideas, the three nodal planes of 
the first layer of mesh I are shown in figures 7, 8, and 9. 

In this work there was only one curved side which was in 
the plane of the triangle, as shown in figure 10. The first, 
seventh, and tenth element nodes were each arbitrarily assigned 
to have as an opposite side the curved side of the triangle. 
The remaining local nodes in each of the respective planes 
of the element were then numbered consecutively in the counter- 
clockwise direction. 

At this point it is appropriate to define connectivity. 

The connectivity of an element is a row vector array that 
relates the local nodes to system nodal points. The connec- 
tivity lists the system nodal points that are "connected" to 
form the element domain in the sequence of local nodal number- 
ing. Thus, the connectivity matrix for mesh I is a matrix 

of size 128 x 15. To illustrate, the connectivity of element 


number 106 is 
S255 5 ClO, 195% 19S, WSs ZIT, 266, lity. bv7,5 ceos We, 10, lily IZ, Soe 


A second finite element mesh (mesh II) consisting of 192 
elements was developed. It contained the same number of ele- 
ments per layer as mesh I. However, mesh II has six layers. 


The second and third layers of mesh I were each divided in 
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Figure 7. Top nodal plane of the first layer of mesh I, 
z = 110 cm 
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Figure 8. Middle nodal plane of the first layer of 
mesh I, z = 95 cm 
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Figure 9. Bottom nodal plane of the first layer of 
mesh I, z = 80 cm 
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Figure 10. Local nodal numbering of curved 
elements 


half, thus forming the two additional layers of mesh II. 
The connectivity matrix and the coordinates of each global 
node for mesh I are given in Appendix A. Appendix B lists 
the connectivity matrix and nodal coordinates of mesh II. 
The reader can construct the different layers and elements 
without much difficulty by using the nodal coordinates and 
the connectivity matrix. A mesh generator was not utilized 


in this work. 


E. COORDINATE TRANSFORMATION 

The use of a quadratic or higher order element permits 
the transformation or mapping of the straight-sided parent 
element into an element with curved sides. Distorted or 
curved elements provide a better fit to curved domains than 
linear elements, and thus a smaller number of elements is 
required to represent the structure adequately. 

The transformation from cartesian coordinates to curvi- 
linear coordinates can be accomplished by employing a one-to- 


one correspondence defined by [11] 


EL 


The element shape functions are utilized to achieve this 


tat Ps 8 | 


1 
2} (46) 
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transformation via the relation 
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where n° is the number of element nodes, and the element 
shape functions N; are in terms of local coordinates. Each 
set of local coordinates corresponds to only one set of car- 
tesian coordinates. 

In performing the transformation, the compatibility re- 
quirement must be met. The transformation into the new, 
curved elements should leave no gaps between adjacent elements. 
If two adjacent elements are generated from parents in which 
the element shape functions satisfy continuity requirements, 
then the curved elements will be contiguous [11]. For the 
isoparametric element, uniqueness of coordinates ensures com- 
patibility. Continuity is assured when adjacent elements 
are given the same sets of coordinates at common nodes. 

Since the element shape functions are in terms of local 
coordinates, an element of volume, dxdydz, must be transformed 
into an element of volume expressed in local coordinates. 

This is achieved through the use of the Jacobian matrix de- 
fined below. Using the chain rule, the relationship between 


Esn,> and a corresponding set of cartesian coordinates 


Seve. iS 


The Jacobian matrix [J] is defined as 
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Jacobian matrix becomes 
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From equation 
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The determinant of the Jacobian is used to transform the vol- 
ume of element in cartesian coordinates to local coordinates. 


For a volume of element [11], 
dxdydz = det [J] d&dndz (51) 


Note that the determinant of the Jacobian is a variable for 
elements of curved geometry. Only in the case of straight- 
sided elements is the determinant of the Jacobian a constant. 


In the plane of the triangle, the area coordinates (Ly, 


L5»b3) number one more than the cartesian coordinates (x,y). 


Thus, L3 is defined as a dependent variable. Tnis establishes 


the origin of the En coordinate system at corner point 3, as 


illustrated in figure 11. Recall that the &n axes 


not be orthogonal. As such 


Using equation (42), 


Applying the chain rule yields 
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Using equations (52a), (52b), and (52c) in equations 
and (53b) gives 


in a triangle 


n& coordinates 
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Now the Jacobian matrix can be evaluated using the shape func- 


tions N. = Ni (L),b,,L4,5) . Inserting equations (54a) and 


(54b) in equation (50) produces 


N,N, 9 aN, aN, ON, 
20L, 4? or U4) , or, : i 
a5 ‘ he oe rao aN hist “ 
J L L wl oS = z = Xs >’ > = y; ’ Dx i Zs 
a3 F w, a3 j ; oly OL; i j oly OL; i 
a Ga ae 
E(——)x. , i(—)y:. , t(——)z. 
9G 1 j la 1 3S a 
(55) 


To summarize, suppose it is required to transform the 


integral 


[ss oF. F'(Lysbosb3.5) dxdydz (56) 
Vol 


to an integral entirely in terms of local coordinates. The 


determinant of equation (55) is then utilized to give 


I oo a F'(Ly sbosbq,)det(d(Ly sb, 5L4,5) Jab dL ode 
=| 0 0 (57) 
Equation (56) is now in a form suitable for numerical inte- 
gration. 
The development of coordinate transformation to this 
point has been under the general assumption that all the 


sides of the element are curved. In this work, the only 
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curved side is in the plane of the triangle. The sides of 
the element along the prism axis are straight. As such, a 
modified Jacobian matrix can be employed, thereby reducing 
the number of calculations to be performed. This modified 


Jacobian is the 2x2 matrix defined by 


oN, N, aN; N. 
oe at a 
> aN, aN, aN, ON; 
[a"] = ate get Se a as 


Along the prism axis or ¢ direction, it can be shown that 


dz = B dz (59) 


where h = height of the element. The volume relationship 


is then given by 


= 


dxdydz = $ det(J"] dljdlys (60) 


The integral of equation (56) then assumes the form of 


} 3 ek 
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Equation (61) is the basis of numerical integration applied 


in this work. 


F. CONSTRUCTION OF ELEMENT MATRICES 


The system matrix operators can be constructed through 


the use of the global basis functions G) or through the 


<n eee 


,s 
_ 


application of the element shape functions. Although the 
global basis functions were used to demonstrate the method 
of Galerkin, the construction of the system matrices in this 
work was achieved through element considerations. The solu- 
tion of the unknown variable w within the element domain 
was approximated by 

+2 kh, ott) 4 te (62) 

ae 

where n° is the number of element nodal points, N. are 
the element shape functions, and v(t) are the time-depen- 
dent nodal magnitudes of ye . The element contribution to 


the system matrix operators is defined by Galerkin's ortho- 


gonality condition expressed by 


J N R™ d¥ob.© ©, $8752. ..0eh” (63) 

Vol 
where R° is defined by replacing b with ye in equations 
(37) and (38), and the integration is over the element volume. 


Using equation (62), the element contribution is portrayed 


for the core by 


e 
CL % NsdVoT] Cp, (t)}-0 ff N,7°N,dVoI]{(c1-y*) 53 uF | N,N Vol ]{(c2+y®) 5} 


Vol 
S ; 
+ LJ NNjen(1 + =~ EN, wy P(t))dvol] {(edeu®) 5} 
td Ylgk Kok f 
+ Cf NjN;avoT] S ett) (cpey(t!) ),dt'} = 0 (64) : 
Vol 


Lar 


ee ew 
> a ae a." 


49 


where dtsdsk = 42,..-505 , and el, ¢2, efc., are constants 
at node i. The bracketed expressions represent square 
matrices of size 15x15, and the braced expressions represent 
column vectors of size 15x1l. For the reflector, the last 
two terms of equation (64) are zero. Before any operation 
can be performed on equation (64), the nonlinear terms (last 
two terms) must be “linearized" and the term with the V2 
operator must be integrated by parts. 15 

The nonlinear feedback term, 2n[l + a ie Nv, (t)] ; 
] 


au : 


was linearized by using predicted values 

of the unknown function at time t. These predicted values 
come from the time integration scheme by Franke [7]. 
Basically, the integration scheme utilizes a predictor-cor- 
rector method which predicts values of the unknown function 
at the next time by using the derivatives of the function. 
Adopting the predicted values a enabled integration over 


space. The term in equation (64) involving the nonlinear 


feedback can be written as 


K: ‘ Pp ; P f P 9 ee } 


The last tern of equation (64) describes the delayed 


neutron contribution. In general, 
t 


c 
oy Th a3 a = es n=] S41 
e Ath f ert we -(t!)dt! =e Mt, ty-1) [e tne] f ." veo(t')dt'] 
0 


=n Mm $e 
+e Ath f ah vai (t!)dt" (65) 
n-1 


where n is number of time steps, t is the current time, 
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and to is the previous time. To approximate the integrals 
in equation (65), the predicted values ve were employed 
in a simple trapezoidal rule. The trapezoidal rule was be- 
lieved to be sufficient since very small time steps were uti- 


lized in this work. For example, at time to 4 


(SUM.) = (previous SUM,)e7*4 + HrerHy (e.)+ yi?) 
j i 2 j ] j 
where 
H = current time step = to - ty 
Xt eee e 
(suM,) =e “"2 f °e veo (t! at! 
0 
ae 
a Ye 
(previous SUM.) = eon x ert w(t! dt! 
0 


The previous sum is formed by accumulating the trapezoidal 


integration of each time step. In general then, for time th 


(SUM,) = (previous SUM. Je” ‘M+ Be My Pct )t¥4) (66) 
The delayed neutron term in equation (64) can be expressed 
as 


[ s NN.dVol] {(¢5+SUM) .} 
Vo 


ee 


In order to bring equation (64) to final form, the V? 


operator was integrated by parts. Since the flux at the sur- 
a4 


face of the reactor is zero, integration by parts yields 


oN. oN, oN; ON. ON, ON, 


2 = a ee SE es RE a A: 
sf N,V N;dxdydz f (x dX is oy oy : oz Zz )dxdydz (67) 
Vol Vol 


In vector notation, equation (67) is expressed as 


aN. 
Bead. 
ox 
: aN. aN aN. aN, 
| NAV N,dxdydz =- f i ame —_ Fra dxdydz (68) 
Vol Vol 
N. 
— 
dz 
Using the chain rule produces 
3x a ox aL, ol 
pe 0 C i aie 
oy oy ay” Ly a 
oN. aN. 
i oc j 
OZ) 0, 0 +3 (=) (69) 


Letting the 3x3 matrix above be (B'], the V* term becomes 


oh My 
( = 
f N,V? Ns dxdydz = aL, aL3 
Vol ON, BN, ON, N,N, ON, aN, 
-f <tr a,)*4ac, * a) seb U8' ITB (ah = ah) paxayaz 
1 2 ws 28 
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where tai}! is the transpose of [B'] . By applying the 
chain rule in equation (47) and using equation (59), it can 


be shown that for this work 


o 
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ca‘! can be derived from equation (71). 
Note from equation (64) that there are three basic ele- 
ment matrices which are defined as follows after applying 


equation (60): 


[655] -3SSS 
00 


* 
NN, detld JdL,dlodz (723 


aN ; 
aa i ] ¥ i oe ay Fk 137 
Sate esd Sr oe a a 
-l10 0 
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Cr om 30, 
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[B'] (5 U5) det[J ]dL,dld¢ 
aN, 
(73) 
wy pak, baby : . A ‘ 
(466, ,] = ISS NAN, 2n( 1+ TN +...4Nyevy6) )detld Jdb dL ode 
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Element matrices [G..] and [6655] are independent of 


ji 
time. However, (6G 55] is time dependent due to the utili- 
zation of the predicted values ¥,° which changes with time. 
In terms of these three basic element matrices, equation (64) 


is 
*@ e e 
[G5 ,]{v;} 1 (GG, ,]icl-v );} + (G,,]i(c2-y );} 


+ [666,,] {(cd-v"),} + [6,,] ((c5+SuM),} = 0 (75) 


iil 
The last two terms of equation (75) are zero for the reflector. 


G. CONSTRUCTION OF THE SYSTEM MATRICES 

The 15x15 coefficient element matrices were calculated 
according to equations (72), (73), and (74); and the results 
were collected element by element into the corresponding sys- 
tem coefficient matrices. The system coefficient matrix 
[BIGG] is developed from (6, ;] , [BIGGG] from [66 55] > 
and [BIGH] from [6665] . [BIGG] and [BIGGG] are inde- 
pendent of time and can be constructed once and for all from 
geometry considerations. [BIGH] is dependent on both geom- 
etry and time due to the time dependence of the predicted 
flux utilized in the feedback term. Thus, [BIGH] is recal- 
culated at each time increment. 

Non-zero contributions to a global nodal point I come 
only from adjacent elements sharing that same nodal point I. 
Thus, the system matrices are sparse and banded. The pro- 


cess of assembling contributions from element matrices 
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requires the identification of a local nodal point (i=1,2, 


..-315) with a global nodal point (I=1,2,...NUMNP, where 
NUMNP is the total number of system nodes). This correspond- i 
ence between element and global nodes is accomplished via 
the connectivity matrix. 
The formal treatment of the field equations in terms of 


the system coefficient matrices is described by the equation 


[BIGG] {vy} + [BIGGG] {(cl-y),} + [BIGG] {(c2-y)/} 


+ [BIGH] {(c4-y),} + [BIGG] {(¢5+SUM),} = 0 (76) 


where the system matrices are NUMNP x NUMNP and the column 
vectors are of length NUMNP x 1. However, the direct appli- 
cation of equation (76) requires a large amount of computer 
storage. To take advantage of the sparsity of the system 
matrices, an optimum compacting scheme described by Ref. 8 
was employed. 

The concept behind OCS is simply to store only the non- 
zero terms of a coefficient matrix. OCS requires two integer 
arrays, say JB and NAME, and a vector of non-zero coeffi- 
cients of the square system matrix. For purposes of illus- 
tration, the square system matrix is called B, and the vector 


of non-zero coefficients of B is called BB. The jth 


integer 
entry in the NUMNP x 1 JB vector is the number q;- This 
number is defined by 

Sy * ik: P; ys VOL yey saws NUNN (77) 


th 


where P; is the number of terms in the i equation. In 
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other words, P; is the number of nodes that the qth node 


“sees". JB is therefore a pointer vector of length NUMNP+1 


whose yeh term locates the initial position in the BB vector 


of the contributing coefficients to the aan equation. The 
NUMNP 
NAME vector of length Mx1l, where M = is P; 
i=] 
NUMNP successive vector blocks of variable length p 


» consists of 


a 
i=1,2,...,NUMNP. The p, integer numbers in the i*” block 

of NAME list the P; contributors to the yeh equation. The 
Mx] BB vector contains the real non-zero coefficients of the 
NUMNPXNUMNP B matrix, arranged in the same contiguous block 
arrangement as the NAME vector. The hak term in the ie 

block or BB(JB(I) + J-1) is B(I,K), where K = NAME(JB(I)+J-1). 
To illustrate, consider the grid shown in figure 12. The two 


array vectors and the coefficient vector matrix of non-zero 


terms are: 
JB = <1, 5, 10, 13, 18, 25, 30, 33, 38, 42> 
NAME = <1,2,4,5/2,3,1,6,5] --- |9,8,6,5> 
BB = <By 1 2B, 9 +BY 4 +By 5! Boo +853 5851 98 az 


267825! ~~~ |Bgg Bg g »Bgg »Bgg” 
In this illustration, NUMNP = 9 and M = 41 [8]. 

In this work, a judicious method of numbering the system 
nodal points was adopted to further reduce computer storage 
requirements. Since the surface boundary nodes of the reac- 
tor represent zero neutron fluxes, the contributions of these 
nodes to interior or non-zero nodes can be discarded. Thus, | 
only the interior nodal points need to be considered. The 


non-zero nodes were numbered first in the finite elemen: 
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Number - element number 


Figure 12. Sample grid used for illustrating OCS 


used so that in the OCS, the number of non-zero nodes (NNZ) 
replaces NUMNP. 

The vectors of non-zero coefficients will be designated 
BIGG, BIGGG and BIGH since the square coefficient matrices 
described in equation (76) were not utilized. To illustrate 
the application of OCS in the system, the following sample 


program is given: 


DO 450 I=1, NNZ 


JBB = JB(I) 
JE = JB(I+1)-1 
DY(I) = 0.0 


DO 500 J=JBB, JE 
LL = NAME(J) 


DY(1) DY(I) + BIGG(J)*H(LL) + BIGGG(J)*c1(LL)* 


W(LL) + BIGG(J)*c2(LL)*p(LL) + BIGH(J)*c4(LL)* 

w(LL) + BIGG(J)*c5(LL)*SUM(LL) (78) 
500 CONTINUE 
450 CONTINUE 


where 
LL = nodal point "seen" by node I 
JE-JBB = total number of nodal points that node I "sees" 
DY(I) = summation of all contributions to node I and 
should sum to zero as stated by the field equa- 
tion 
The assumption of a homogeneous reflector was relaxed in the 
application of equation (78). Interface nodes were assigned 


properties of the core. Therefore, if LL is a reflector 
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node not on the core reflector interface, the last two terms 


of equation (78) are nonexistent. 
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TV. NUMERICAL INTEGRATION 


A. LINE AND AREA INTEGRATION 

The solution of the element matrix equations was achieved 
through numerical integration since an exact closed form solu- 
tion cannot be established. The volume integration was accom- 
plished by using a line integration in the g-direction and an 
area integration in the plane of the triangle. The line inte- 


gral is described by the Gaussian quadrature formula [12] 
] 


n 
Jf eteyac = 2, a fla,) (79) 
=) k=] 


where n is the number of Gauss integration points 


= 
iT] 


weighting coefficients 


the function f(c5) evaluated at Gauss point ap 


Table III lists ays Hy, and n ek oe 


k 
The area integration was achieved by the equation 


1 ,leL 
] m 
~ m m m 
ry f(L,,L,,L,)dL,dL, = 2 waf(L, sly Lb, ) (80) 
0 ie} 


m=] 


where m is the number of area integration points and Wa 
are the weights. The numerical integration points for the 


area integration are given in Table IV which was extracted 


TABLE III 


Abscissae and Weight Coefficients 
of the Gaussian Quadrature Formula 


] n 
I f(c)dc = DY Hy F(a,) 
=f k=1 
ta H 
n=2 
057735 02691 89626 100000 00000 00000 
: n=3 
0:77459 66692 41483 0:55555 55555 55556 
0:00000 00000 00000 0:88888 88888 88889 
n=4 
086113 63115 94053 034785 48451 37454 
0:33998 10435 84856 065214 51548 62546 
n=5 
090617 98459 38664 0:23692 68850 56189 
0:53846 93101 05683 0:47862 86704 99366 
000000 30000 00000 0:56888 88888 88889 
n=6 
0-93246 95142 03152 0-17132 44923 79170 ’ 
0-66120 93864 66265 0:36076 15730 48139 
0:23861 91860 83197 0:46791 39345 72691 
n=7 
094910 79123 42759 0-12948 49661 68870 
0:74153 11855 99394 0:27970 53914 89277 
040584 51513 77397 038183 00505 05119 
0-00000 00000 00000 041795 91836 73469. 
n=8 
0:96028 98564 97536 0-10122 85362 90376 
0:79666 64774 13627 0:22238 10344 53374 
0:52553 24099 16329 0:3!370 66458 77887 
9:18343 46424 95650 0:36268 37833 78362 
n=9 
0:96816 02395 07626 008127 43883 61574 
083603 11073 26636 018064 81606 ©4857 
0-61337 14327 00590 9-26061 06964 u2935 
032425 34234 03809 031234 70770 40002 
0:00000 00000 00000 033023 93550 0!260 
n= 10 
097390 65285 17172 0:06667 13443 08688 
086506 33666 88985 0:14945 13491 50581 
067940 95682 99024 0:21908 63625 15982 
043339 $3941 29247 0:26926 67193 09996 


14887 43389 81631 0:29552 42247 14783 


TABLE IV 


Numerical Formulas for Triangles 


Order Error Points Tnangular Weights 
Co-ordinates 2W, 

Linear R = (h*) a 4.4.4 1 

a 44,0 4 
Quadratic R = Oh) b 0.4.4 | 

e «6. $, 0,5 } 

a 4.4.4 -h 
Cubic R = Mh*) 6 HAA Hi 

co HL 

d *&. AH 


ae 2 
bh 4ho 
¢ 0.4.4 
R = Oh‘) d $0.4 
e 1.0.0 
t 0.1.0 
g 0.0.1 


| Cubic : 

; 

‘ Pe ee ne ae ee ee 
{ 


a hg 0-225 

b x1. By. B, 

¢ By. 2,. B, 013239415 
Quintic R = Kh‘) d Brnbys 

e 1;. By, B; 

f By. 23. By 0:12592918 

9 as Zp 


Bp Ppp. 


with 

2, = 005971587 
B, = 047014206 
2, = 079742699 
B, = O1012R6S1 
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from Ref. 11. The volume integration of the function 


F(L,,L,5b350) using numerical integration is therefore 


11 1-Ly 
-“F 6 0 
n m 
& Wy Ze Wf (eyes t gM sk gts) (81) 
= m= 


where Mec Hy . In the manner of equation (81), the element 


matrices become 


n m 
Pee S, Ww LK cr" 
(65 ,] a3 my We x Wo Nj (Ly rLy rb, ot )N,det{d ] (82) 
rec.) = = > wow F (Ly LM Lat 08) det(J] (83) 
Fe ee yy eg 
n m * 
eee <5 ay Sw, Thy ey tye) deste | 484) 


ro 
it 
at 
| 
“ 
pe 


where F -and T are easily derived from equations (73) and 
(74), respectively. Note that N. and det[J ] are also 
evaluated at each integration point. They are not shown as 
such merely for the sake of convenience in writing the equa- 


tions. 


B. NUMBER OF INTEGRATION POINTS 

It is difficult to estimate the number of integration 
points required for good accuracy due to the complexity of 
the functions involved. The basic rule that the best number 
of integration points is found by trial and experience was 


adopted. For the (65;] element matrix, the numbers involved 
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can be approximated. This was done by letting the determi- 
nant of the Jacobian equal to twice the area of the curved 
triangle (checking det[J*] at each integration point showed 
that this assumption was not too unreasonable). With det[J*] 
outside the integration process, (65,] can be solved in 

closed form by integrating out ¢ from -1 to +1 and then apply- 


ing the closed form equation [12] 


m m m m,!m,!m.,! 
] 2 3 - EPEC 
3 d(Area) = 2A m, #m,#m4¥2 


Area 
(85) 
where My, M5, Mz are positive integer exponents and A _ is 
the area of the triangle. 
Five test points within the 15x15 element matrix were 
selected, as listed in Table V. The values obtained from 
the application of equation (85) are also given in Table V. 
Three sets of area integration points, each with a different 
number of c¢ Gauss points, were used. These three sets of 
area integration points, given in Table IV, are: 
1. cubic order (4 points) 
2. cubic order (7 points) 
3. quintic order (7 aortas) 
Using each of the three integration points above with differ- 
ent c Gauss points in equation (82) produced the results 
obtained in Table V. From these results, the quintic order 
area integration points were selected for the element matrix 


(6, ;] with three <¢ Gauss points in the prism axes. This 


set of integration points was also used for [GGG 5]. 


eee ee 


ae eae Nt ae 
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TABLE V. 


Selection of Integration Points for [G55] 


Area 

cS points oints G(1,1 G(2,2 G(1,9 G(6,9 G(9,9 

13 cubic 1415 5278 -2756 5072 8536 
(4 pts) 

5 us 1817 5278 -3170 5072 10240 

7 L 1817 5278 -3170 5072 10240 

3 cubic 3248 8247 -3114 4960 10243 
(7 pts) 

5 " 3249 8247 -3114 4960 10240 

7 " 3249 8247 -3114 4960 10240 

3 quintic 2464 6600 -3144 5020 10240 
(7 pts) 

5 * 2464 6600 -3144 5020 10240 

7 Ms 2464 6600 -3144 5020 10240 

Approximated values 2500 6700 -3142 5026 10053 


~~ 


The numbers for the (6G, 5] element matrix could not 
be approximated. The best that could be done was to obtain 
an idea of the order of magnitude of this element matrix. 
Towards this end, a linear triangular element was assumed. 
Using linear approximation, the order of magnitude was found 
to be about 109, Using the three different area integration 
points mentioned above with varying ¢ Gauss points in equa- 
tion (83) yielded the results given in Table VI. Further 
checks on the [6G 55] element matrix showed that the cubic 
order with four area integration points yielded the desired 
order of magnitude. Thus, the fourth order cubic with five 
5 Gauss points was employed for the (GG. ,] element matrix. 
A note should be mentioned here in regards to the vast dif- 
ference in results obtained for (4G, ,] using different 
area integration points. Most likely, it was due to the 
[B'] and [B'y! matrices which required the inversion of 
sr, ’ a, » etc. Or perhaps it was caused by the nature of 


the hybrid element used in this work. In any case, further 


investigation is warranted in this area. 


ct points 
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TABLE VI 


Selection of Integration Points for (665, J 


Area 
oints 


cubic 
(4 pts) 


cubic 
(7 pts) 


quintic 
(7 pts) 


GG(1,9 
18.4 


22.0 
22.0 


-1,03x10* 2.1110 


-1.03x10" 


-1.03x10" 


-52.6 


-52.6 


-52.6 


6G(3,3) GG(14,9) GG(8,8 
170.1  -165 165.9 
177.6  <="86 195.9 
177.6. -186 195.9 

7 1.05x10” 8.42x107 
2.11x10” 1.05x10’ 8.42x107 
2.11x10’ 1,05x10’ 8.41x10/ 


1.27x10° -6.64x10"° 6.84x10° 


1,27x10* -6.64x10* 6.84x10° 


4 4 


6.84x10° 


1.27x10° -6.64x10 


GG(15,15 
106.1 


106.1 
106.1 
2056 


2056 


2056 
1.51x10° 
1.51x19° 


1.51x10° 


V. TEST PROBLEMS AND RESULTS 


tions in the form of a ramp input described by 
* 
Ue = Xe +. at (86) 


| The reactor was subjected to uniform and local perturba- 


where a is the change in re per unit time and is 
the critical fission cross-section. arg must first be ob- 
tained before the perturbations can be applied. This was ac- 
complished by trial and error until a stationary solution 

was reached. For mesh I, I, was found to be 0.0057360 

per cm. No attempt was made to find ay for mesh [I due 

to time limitations. As such, the test problems outlined be- 


low were applied to mesh I. Future work is planned to apply 


En 


the test problems on mesh II. 


The following perturbations were applied: 


a) Uniform perturbation of 10 dollar of reactivity per 


second: 
Ze(r,t) = leg = Ot in the core 
where a = 0.005893/cm-sec 


b) Local perturbation at the core center of 100 dollar 


of reactivity per second: 


* _~ : 
i,(r,t) ae? +atd(r,) » in the core 


of 
—_— , 68 j 


t 
| 


See es one 


where . is (0,0,0) and a = 0.015407/cm-sec 


c) Local, off-center perturbation: 


Z,(r,t) Se a.td(r,) ; + = 1,253, in the core 


where 
ry = (0,60,40) 
a, = 0.015407/cm-sec = 100 dollar per second 
do = 0.008123/cm-sec = 50 dollar per second 
Ge = 0.005894/cm-sec = 10 dollar per second 


Three test points, (0,0,0), (60,0,0), and (-60,0,80), 
were selected to trace the neutron time history. For cases 
a) and b), the neutron flux was plotted at each test point 
during transience. This is shown in figures 13 and 14. 

Case c) involved three ramp inputs and were conducted for 
both the linear and nonlinear reactor equations. The linear 


and nonlinear responses were compared at each test point for 


each ramp input and are illustrated in figures 15 through 23. 


The radial and axial flux distributions at time t = 0.0123 
second were plotted for the steady state and the 100 dollar 
perturbation case. These are embodied in figures 24 and 25. 
Finally, a@ neutron flux early time history between mesh I 
and mesh II was plotted to check the effect of using a finer 
element mesh. The result is portrayed in figure 26. 

Figures 13 and 14 revealed a clear space dependence of 
the neutron flux during transience as expected. Figures 15 


through 23 demonstrated the effects of temperature feedback 


Se ee eee 
> ‘ot ch 


only the interior nodal points need to be considered. These 


non-zero nodes were numbered first in the finite element mest 
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A - neutron flux at test point (0,0,0) 
B - neutron flux at test point (60,0,0) 
C - neutron flux at test point (-60,0,80) 


Neutron flux transient history at various 
test points for a local central perturbation 
of 100 dollar/sec of reactivity. 
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Figure 15. Linear and nonlinear fluxes at (0,0,0) due 
to a local perturbation of 100 dollar/sec of 
reactivity at (0,60,40). 
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Linear and nonlinear fluxes at (60,0,0) due 
to a local perturbation of 100 dollar/sec at 
(0,60,40). 
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Figure 17. 
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A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (-60,0,80) 
due to a local perturbation of 100 dollar/sec 
at (0,60,40). 
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Figure 18. Linear and nonlinear fluxes at (0,0,0) due 
to a 50 dollar/sec local perturbation at 
(0,60,40). 
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Figure 19. Linear and nonlinear fluxes at (60,0,0) 
due to a 50 dollar/sec local perturbation 
at (0,60,40). 
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Figure 20. Linear and nonlinear fluxes at (-60,0,80) 
due to a 50 dollar/sec local perturbation 
at (0,60,40). 
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Figure 21. Linear and nonlinear fluxes at (0,0,0) due 
to a 10 dollar/sec local perturbation at : 
(0,60,40). ; 
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Figure 22. Linear and nonlinear fluxes at (60,0,0) due 
to a 10 dollar/sec local perturbation at 
(0,60,40). 
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Figure 23. Linear and nonlinear fluxes at (-60,0,80) 
due to a 10 dollar/sec local perturbation 
at (0,60,40). 
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Figure 24, Radial flux distribution for the steady state 
and 100 dollar/sec local perturbation. 
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Figure 25. Axial flux distribution for the steady state and = 
100 dollar/sec local perturbation. a 
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Figure 26. Early time history of the neutron flux at core 
center using mesh I and mesh II. 
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and delayed neutrons on the flux. To obtain an idea of the 
difference in magnitude between the linear and nonlinear 

flux for the 100 dollar local perturbation, at time t = 10 
second, the linear case predicted a neutron flux at the core 


center of 2.024 x 10!2 neutrons per cm? per sec. The non- 


linear reactor predicted 1.85 x 19! neutrons per em? per 
sec. Although the numbers are not sufficiently accurate due 
to the crudeness of the finite element mesh, the significant 
difference in the order of magnitude between the linear and 
nonlinear neutron flux leads to the belief that the three- 
dimensional finite element model utilized is predicting the 
correct trend of neutron flux behavior. 

The radial flux distribution for the steady state shown 
in figure 24 portrays the expected flux distribution. The 
radial flux distributions for the local perturbation case 
show a skew distribution at the point of perturbation. The 
axial flux distribution at steady sta e is very much symme- 
tric about the center. Again, the expected skew at the 
point of perturbation is there, as shown in figure 25. 


* 
Figure 26 gives an indication that the © for finer 


mae 
meshes is higher. However, curve B of figure 26 should be 
extended to longer times to verify this hypothesis. Curve B, 
as plotted in figure 26, used four hours of computer time 


employing the H-compiler of the IBM 360/67. 


VI. CONCLUSIONS AND RECOMMENDATIONS 


The finite element mesh employed here is crude, and, 
thus, results that were obtained should be considered as posi- 
tive indicators rather than numerically conclusive facts. The 
trend of neutron flux behavior is the major thrust of this 
work. From the results obtained, it is concluded that the 
expected patterns of neutron behavior as predicted by the 
three-dimensional finite element model used do occur. These 
patterns were best demonstrated by the differences between 
the linear and nonlinear flux responses and by the spatial 
flux distribution at steady state and during local perturba- 
tion. The three-dimensional quadratic finite element model 
utilized should produce better results by resorting to a 
finer mesh. The draw~back to a finer mesh, of course, is 
the significant increase in computer time and storage re- 
quirements. Once accurate results are obtained through 
finer element meshes, comparisons between three- and two- 
dimensional models can be attempted. 

A mesh generator was not developed for this problem. It 
is recommended that this be done to ease the transition from 
one mesh to another and to minimize human error. In addi- 
tion, a similar calculation using a three-dimensional linear 
element should be performed to corroborate the results ob- 
tained here. It should be particularly noted that this type 


of problem is highly sensitive to the fission cross-section. 
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In the search for ae » it was necessary to adjust le to 
the sixth decimal place. There is no exact method of deriv- 
ing ce due to the highly nonlinear aspect of the problem; 
therefore, trial and error must be used. Finally, the 
Gaussian quadrature used to determine the element matrices 
should be investigated. The cause of the differences in 
values obtained by using different number of integration 
points should be established. Was it due to the integration 
points, the coordinate transformation, or the element shape 
functions? This is an important question upon which future 


numerical results could be based. 
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APPENDIX A 


MESH I CONNECTIVITY AND COORDINATES 


The connectivity matrix for the 128-element mesh is as follows 


11 12 13 14 15 
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